function uh = SiWeno_Limiter(uh)

global Nx dimPk bcL bcR NumGLP phiG weight mm

uhb = [[0,0,0];uh;[0,0,0]];
uhmod = zeros(Nx,dimPk);
uhmod(:,1) = uh(:,1);
temp=zeros(Nx,NumGLP);
uhG = zeros(Nx,NumGLP);
tol=10e-6;
for i = 1:Nx
    for d = 1:dimPk
        uhG(i,:) = uhG(i,:) + uh(i,d)*phiG(:,d)';
    end
end
uhbx = [[0,0,0];[0,0,0];uh;[0,0,0];[0,0,0]];
uhGb=[[0,0,0,0,0];uhG;[0,0,0,0,0]];

% set_bc
if bcL == 1
    uhb(1,:) = uh(end,:);
    uhbx(1,:) = uh(end-1,:);
    uhbx(2,:) = uh(end,:);
    uhGb(1,:)=uhG(end,:);
end

if bcR == 1
    uhb(end,:) = uh(1,:);
     uhbx(end,:) = uh(2,:);
     uhbx(end-1,:) = uh(1,:);
    uhGb(end,:)=uhG(1,:);
   
end
 
for i = 1:Nx
    deltaUR = uh(i,2) + (2/3)*uh(i,3);
    deltaUL = uh(i,2) - (2/3)*uh(i,3);
    deltaURM = uhb(i + 2,1) - uhb(i + 1,1);
    deltaULM = uhb(i + 1,1) - uhb(i,1);
    
    deltaURM1 = minmod(deltaUR,deltaURM,deltaULM);
    deltaULM1 = minmod(deltaUL,deltaURM,deltaULM);
  
    uhmod(i,2) = (deltaURM1 + deltaULM1)/2;
    uhmod(i,3) = 3*(deltaURM1 - deltaULM1)/4;
end
%uh=uhmod;
for i=1:Nx
   if(uh(i,2)-uhmod(i,2)>tol||uh(i,3)-uhmod(i,3)>tol)
       
      for j=1:NumGLP
          temp(i,j) = WENO5(uhGb(i,j), uhGb(i+1,j), uhGb(i+2,j), ...
                  uhbx(i,1), uhbx(i+1,1), uhbx(i+2,1), ...
                  uhbx(i+3,1), uhbx(i+4,1));

      end 
     
        uh(i,:)=0;
        for d = 1:dimPk
             for i1 = 1:NumGLP
                 uh(i,d) = uh(i,d) + 0.5*weight(i1)*temp(i,i1)*phiG(i1,d);
             end
        end

            for d = 1:dimPk
              uh(i,d) = uh(i,d)/mm(d);
            end
    end
end

 end

function a1 = minmod(a,b,c)

global hx M

if abs(a) < M*hx^2
    a1 = a;
else
    if sign(a) == sign(b) && sign(a) == sign(c)
        a1 = sign(a)*min(abs([a,b,c]));
    else
        a1 = 0;
    end
end

end

function h = WENO5(a,b,c,k1,k2,k3,k4,k5)

epsilon = 1e-6;

h0 = a-k2+k3;
h1 = b;
h2 = c-k4+k3;

beta0 = 13/12*(k1 - 2*k2 + k3)^2 + 1/4*(k1 - 4*k2 + 3*k3)^2;
beta1 = 13/12*(k2 - 2*k3 + k4)^2 + 1/4*(k2 - k4)^2;
beta2 = 13/12*(k3 - 2*k4 + k5)^2 + 1/4*(3*k3 - 4*k4 + k5)^2;

w0 = 0.01/(epsilon + beta0)^2;
w1 = 0.98/(epsilon + beta1)^2;
w2 = 0.01/(epsilon + beta2)^2;

S = w0 + w1 + w2;
w0 = w0/S;
w1 = w1/S;
w2 = w2/S;

h = w0*h0 + w1*h1 + w2*h2;

end
